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1 Introduction 

In this paper we analyze geometric numerical integration algorithms for weakly 
dissipative perturbations of Hamiltonian systems. The dissipation rate is gov- 
erned by a small parameter < e ^ 1, as in the simple example 

^ = V (1) 

p ^ -V {q) - eDq ^ ' 

where q,p G R", M is a symmetric positive definite matrix, and D is a sym- 
metric positive semi-definite matrix. At e = the system is Hamiltonian with 
H{q,p) = -^p^ AI^^p + V{q), and for e > it is dissipative. This type of 
dissipation is the simplest example of Rayleigh damping (the general case is 
presented in Section |2] below) . 

There are many applications. For example, in rolling bearing design, the ob- 
jective is to minimize friction losses. The dynamics of such mechanical systems 
is often well described by this model. In numerical simulations, it is then of im- 
portance that energy losses are correct. This means that a geometric integrator 
(nearly) conserving energy at e = 0, and giving the (nearly) correct dissipation 
rate as a function of e > 0, is desirable. Moreover, if possible, explicit meth- 
ods are preferred. However, well-known reversible explicit methods, such as the 
two-step explicit midpoint rule, typically lose their stability in the presence of 
dissipation. Furthermore, explicit symplectic methods for separable Hamilto- 
nian systems (e.g. Lobatto IIIA-IIIB) become implicit if dissipation is added. 
This motivates a special study and analysis of explicit geometric integrators for 
systems with weak Rayleigh damping. 

The paper will address this question by examining splitting methods that 
take advantage of the linearity of the dissipative perturbation. Three integra- 
tion schemes are investigated. All of these schemes are explicit in the sense that 
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no numerical root solving is required. The first two of the methods are expo- 
nential integrators (meaning they require the solution of a linear ODE). The 
third method is either fully explicit, or linearly implicit (meaning it requires the 
solution of a linear equation). In our numerical examples, the three methods 
under study turn out to yield practically identical results. Although the pre- 
sented methods are all second order accurate, it is an easy extension to construct 
higher order methods by using higher order compositions. By backward error 
analysis we show that the methods have asymptotically correct dissipation be- 
havior as the dissipation parameter e — ^ 0"*". Moreover, we analyze the decay of 
the modified Hamiltonian function, and give theoretical results on monotonicity 
of its evolution for small enough step sizes. We also analyze conservation of 
momentum, and near conservation of relative equilibria, in the case when both 
the conservative and the dissipative parts are invariant under a symmetry. 

The examined methods are invariant under choice of coordinates in the con- 
figuration space. Therefore, they are well defined on any smooth configuration 
manifold without reference to a specific choice of coordinates. Throughout this 
paper we utilize the differential geometric framework of Riemannian geometry to 
describe the systems of interest. Although an important consequence, the main 
reason for doing so is not to include systems on non-linear spaces (constrained 
systems). Rather, the presented framework and corresponding analysis become 
more transparent in a differential geometric setting. Also, the fact that the re- 
sult are derived for abstract manifolds is essential for some of the analysis carried 
out in the paper. (For example, independence of configuration coordinates is a 
key ingredient in our proof of preservation of symmetry invariance.) Moreover, 
well known results in differential geometry, in particular on structural stability, 
symmetry groups and momentum maps, and relative equilibria, become more 
easily available in this setting. For the convenience of the reader, we sometimes 
carry out invariant calculations even in cases when they are evident from a 
differential geometric point of view. 

The idea of applying geometric integration schemes to weakly dissipative 
problems is certainly not new. For instance, it is well known that symplectic 
methods applied to perturbed integrable systems (e.g. the Van der Pols equa- 
tion) are superior to standard methods, in the sense that weakly attractive 
invariant tori are much better preserved; see |H1 Chap. XII]. (Although the con- 
cepts of relative equilibria and invariant tori are closely related, the latter can be 
seen as a special case of the former, the settings and analyses in this paper and 
in [HI Chap. XII] are not overlapping. A comparison between the two settings 
is given in Section |6.1| below.) Furthermore, in the paper |12| it is observed 
numerically that variational methods, in particular a semi-explicit version of 
the Newmark method, applied to weakly dissipative systems give much more 
accurate energy dissipation rates than standard methods. In the simplified case 
of conformal Hamiltonian systems, corresponding to _D = Id in equation Q 
above, it is possible to construct integrators that exactly preserve the corre- 
sponding geometric structure; see |18l [20]. Concerning preservation of relative 
equilibria and numerical integration, there has been extensive work within the 
framework of energy-momentum methods; see [S1IH[2|3]. Under suitable condi- 
tions, energy-momentum methods exactly preserve relative equilibria. However, 
energy momentum schemes are fully implicit. Also related to our work is the 
recent paper |14| . which uses a splitting technique, similar to the one used in 
this paper, for the integration of post-Newtonian equations, possibly with weak 
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dissipation due to relativistic effects. 

The paper is organized as follows. A problem description is given in Sec- 
tion [2] Numerical integration schemes are presented in Section [3] A backward 
error analysis of the methods is given in Section |4] The result of the analysis 
is verified by a numerical example (planar pendulum) in Section |4.3[ In Sec- 
tion [5] we show that the methods preserve symmetry invariance, and conserve 
corresponding momentum maps. In conjunction, Section [6] contains a study of 
conditions for near conservation of relative equilibria, verified by a numerical 



example (elastic pendulum) in Section 6.2 Finally, conclusions are given in 
Section [T] 

Throughout the paper we assume that the reader is familiar with basic dif- 
ferential geometry. In particular, the following geometric concepts are used 
without reference: Riemannian manifolds, tangent and co-tangent bundles, ten- 
sor fields. Lie derivative, symplectic maps, Hamiltonian vector fields, flow of a 
vector field, exponential map. For the convenience of the reader, other more 
advanced, but still fundamental, geometric concepts are introduced and defined 
as they are used. The introductory textbook on geometric mechanics by Mars- 
den and Ratiu [T7] is more than enough to cover all the differential geometric 
concepts in this paper. 

We adopt the following notation, i? denotes a manifold, and T*£2 its 
corresponding tangent and co-tangent bundle. Throughout the paper q denotes 
an element in =S. Furthermore, p denotes an element in the co-tangent space 
T*^. Thus, p is always associated with q. The pairing between an element p S 
T*^ and g G is denoted (p, q). We use z for a point in T*^ corresponding 
to {q,p). The derivative of a function (/) on ^ at the point q (also called the 
tangent map) is denoted Tqcj). The set of smooth vector field on a manifold ^ 
is denoted and the set of diffeomorphisms of ^ is denoted Diff(^). 

The flow of a vector field X G is exp{tX), where t is the time parameter 

and exp : X{^) — >■ Diff( J^) is the exponential map. The Lie derivate along X is 
denoted £x- The space of smooth real- valued functions on ^ is denoted 3^{£P). 
If S 5'(T*^), then its corresponding canonically Hamiltonian vector field 
is denoted Xjj. The canonical Poisson bracket between H,G G 3^{T*£2) is 
denoted {H, G}. The space of type (r, s) tensor fields over ^ is denoted (^). 



2 Problem Description 

Let (i?, g) be a Riemannian manifold of dimension n. where g e 7^(^) is 
a metric tensor, i.e., symmetric and positive definite. In our context ^ is 
the configuration space of a conservative mechanical system with Lagrangian 
function of standard form L{q,q) — ^gq{q,q) — V{q). Thus, the metric tensor 
determines the kinetic energy of the system. 

We are interested in dissipative perturbations of such systems. To this end, 
let d € 'J^(cS) be a symmetric positive semi-definite tensor field called dissi- 
pation tensor. We shall assume that the rank of d is constant, i.e., that it 
is a regular tensor. Further, the Rayleigh dissipation function is defined by 
R(q,q) — 5dq(g,<7). We consider the following dissipative perturbation of the 
Euler-Lagrange equation 

ddL _dL _ ^^dR 
dt dq dq dq ' 
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where e > is the perturbation parameter. 

From a structural, as well as a numerical, point of view it is often prefer- 
able to transform the Euler-Lagrange equation into a first order Hamiltonian 
system on the co-tangent bundle phase space T*^. The metric tensor induces 
an isomorphism M(g) : T^J? — )■ T*i2 by the flat operator c^' := gq{q, •). Its 
inverse defines the sharp operator p i— > :— M.{q)^^p. M is sometimes called 
the inertia operator. Since the kinetic energy descends from the metric tensor, 
the Legendre transformation is given by the fiat operator, so the momentum 
variable is p :— dL/dq — q^ . The Hamiltonian function, describing the total 
energy of the system, is given by the sum of kinetic and potential energies 

Hiq,p) = lgg{pKp^)+Viq). (3) 




Substituting variables, the perturbed Euler-Lagrange equation ^ trans- 
forms into a perturbed Hamiltonian system 

P=-'^-edg{p^,-), ^ ^ 

where we have used that dH/dp = p", which follows since gq{p^,p^) = {p,p^)- 

Since the Hamiltonian part of equation ^ is perturbed by a non-Hamiltonian 
perturbation, total energy is not conserved. Indeed, the time evolution of the 
total energy along solution curves is given by 

d , dH . dH 



-H{q,p) = -^q+^P = -£d,(p«,p«) = -edq{q,q) < 0. (5) 



Thus, the role of the dissipation tensor is to specify the rate of energy decay at 
each point in phase space. 

The form of system Q is invariant under point transformations, as we will 
now derive. Indeed, let ^ and be diffeomorphic manifolds and </> : ^ — ^ =S a 
diffeomorphism. Consider now the change of variables (s, r) — T*(j){q,p), where 
T*(j) : T*^ T*y is the co-tangent lift of 4) (cf. Marsden and Ratiu [T71 
Sect. 6.3]). (In coordinates, (s,r) o is defined by the generating function 

of second kind G-ziq, r) = (j){q)-r.) We know already from standard Hamiltonian 
theory that the governing equation for the Hamiltonian part of the vector field 
in the variables (s,r) transforms into a new canonical Hamiltonian system for 
a new Hamiltonian function K on T* given by 

K{s,r)^H{T*4-Hs,r)) 

= (r, (r,^)-iM(0(s))-i(r,0)-i*r) + 

V ' 

= (r,r«) + (0*y)(s) 

= (rg).(r«,r») + (0*F)(,s) 

where {(j)*g)s{u,v) = g^(s){Ts(j)u,Ts(l)v) and {(l)*V){s) = V{4>{s)) is the pull-hack 
of g and V . Thus, in the new variables (s, r) the Hamiltonian part of system Q 
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takes again the same form but with the pull-backed quantities 0*g and (l)*V. 
Next, the vector field over T*^ corresponding to the dissipation is given by 
eY{q,p) = (0,-ed,(p«,-)) (0, -eD(gy), where : r,i2 T*^. The 

corresponding vector field in the (s, r) variables is the push-forward of Y by r*(/), 
which is given by (r*0)*y := T{T*4>) oXoT*(j)-^. In matrix notation 



(r»,r(s,r) = 



-1 





B(0(s))(r,(/.)-i*r 



(r,^)*2)(0(s))M(^(s))-i(T,0)-i*r 



(6) 



(r,0)*D(0(5))r,0(r,(/))-iM(0(s))-i(T30)-iv 
V V ' 

(0*d),(..) r« 

= (O,(0*d).(r«,.)) 

In summary, this means that system Q is determined by the three geometric 
objects (tensors) g, d and V in such a way that a change of variables (s, r) = 
T*(j)(q,p) takes system ^ into a new system of the same form defined by the 
pulled back triple 0*g, (f>*d and (j)*V. In other words, the following diagram 
commutes: 



g,d,V 



fiow 



q{t),p{t) 



h,e, 



flow 



s{t),r{t) 



(7) 



Thus, the form of equation Q is independent of the choice of coordinates 
on The form is not, however, invariant with respect to any symplectic coor- 
dinate transformation on the full phase space T*cS: only point transformations 
are allowed. As we will see, independence of coordinates on =S is preserved by 
the numerical schemes suggested in Section [3] 



3 Numerical Integration Schemes 

In this section we present the explicit and semi-explicit time-stepping schemes 
for the system Q above to be analyzed in the paper. To begin with, we rewrite 
equation Q in the more compact form 

z = XT{z)+Xv{z)+eY{z), z^{q,p)eT*^, (8) 

where Xt^Xy are the Hamiltonian vector fields corresponding to T and V 
respectively, and eY is the non-Hamiltonian Rayleigh dissipation part. 

Since the potential energy function V is independent of p, the fiow exp(<Xy) 
of Xy is explicitly integrable for any choice of coordinates on ^ simply by ap- 
plying the forward Euler method. Throughout the rest of this paper we make 
one assumption which is crucial for the implementation of the suggested meth- 
ods: that we know coordinates on ^ in which the fiow expliXx) is explicitly 
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computable. For particle systems, this is typically accomplished by choosing 
Cartesian coordinates in which the inertia operator M is independent of q. In 
this case, exp(tXj-) is again computable using the forward Euler method. 
We study the following schemes: 

Algorithm 3.1 (Three-term splitting method). Let /i > be the step size and 
define ^ z^+i by the map 

<I>(jg ^ cxp(^£y) o exp(^XT) o eyip{hXv) o exp(^XT) o exp(^irr). 

Algorithm 3.2 (Two-term splitting method). Let /i > be the step size and 
define Zk ^ Zk+i by the map 

= exp(^XT) o exp {h{Xv + eY)) o exp(^XT). 

Algorithm 3.3 (Runge-Kutta splitting method). Let h > he the step size 
and let be a Runge-Kutta method for X e X(r*=S). Define Zk t-^- Zk+i by 
the map 

*RKS = exp(^XT) o o exp(^XT). 

All of these methods reduce to the classical Stormer-Verlet method $gy 
(also known as the leap-frog method or the Verlet scheme) in the case e = 0. 
(To be more precise, they reduce to the dual version, or _B-version, of the clas- 
sical Stormer-Verlet method; see [7 for various interpretations of this method.) 
Furthermore, Algorithm |3.1| and Algorithm |3.2| are second order accurate, as 
is Algorithm |3.3| if the Runge-Kutta method used is at least of order 2. It is 
straightforward to extend these methods to higher order, by using higher order 
compositions (see |19| for a review of splitting methods). 

It should be noted that the vector field Xy + eY is linear in p and thus 
explicitly integrable using the exponential Euler method (cf. [10 ). Hence, 
explicitly computable. Likewise, Y is explicitly integrable using the exponential 
Euler method , so is also explicitly computable. If the Runge-Kutta method 



in Algorithm 3.3 is explicit, then ^p^j^S ^® f^^^y explicit. If the Runge-Kutta 
method used is implicit, then each step of •I'j^j^g requires the solution of a linear 
system. 

Runge-Kutta methods are well defined on linear spaces. Thus, one might 
think that Algorithm 3.3 is defined only in the case when T*J2 ~ M^". This is 



however not the case, since the dynamics of Xy + sY is trivial in the configu- 
ration variable q, so each "Runge-Kutta step" '^Xv+sY takes place entirely on a 
co-tangent space T* £2, which is a linear space. Furthermore, since Runge-Kutta 
methods are invariant with respect to linear coordinate changes, and since any 
smooth change of configuration coordinates corresponds to a linear change of 
coordinates on each co-tangent space. Algorithm [33] is defined intrinsically (i.e., 
without reference to a specific choice of configuration coordinates). 

Remark 3.1. Notice that system Q is of the form 

q = f{q,p) 
p = 9{q,p) 
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so a partitioned Runge-Kutta method may be used to integrate it. The classical 
Stormer-Verlet method can be extended to the partitioned Runge-Kutta method 
given by the second order Lobatto IIIA-IIIB pair; see Sect. II. 2]. From 
our point of view, assuming that we are using a coordinate system in which 
the kinetic energy is independent of q, this extension of the classical Stormer- 
Verlet method is conjugate to Algorithm |3.3| with the implicit midpoint rule as 
choice of Runge-Kutta method. Indeed, if and $g denotes, respectively, 

the implicit midpoint rule, the implicit Euler method a nd t he explicit Euler 
method, it holds that = ^'l^^ o Thus, Algorithm [sisl with the impHcit 



olds that <I>{^ = *E o *i • Thus, Algorithm 

midpoint rule as the choice of Runge-Kutta method is exp^hXr /2)o^'^^^ o^'^^^ o 

exp(ft.XT/2). This method is conjugate to ^J*^^ o exp(/iXT) o <i>g^'^, which is 
exactly the A-version of the partitioned Runge-Kutta extension of the classical 
Stormer-Verlet method. 



As mentioned above, the suggested integrators are invariant with respect to 
choice of coordinates on the configuration space. From a differential geometric 
point of view, this result is obvious, since all the objects and operations used 
in their definition are intrinsic (configuration coordinate independent). How- 
ever, for the convenience of the reader, we carry out a proof with complete 
calculations. 



Proposition 3.1. Algorithms \3.1W3.3\ are invariant with respect to choice of 
coordinates in £S. That is, if 12 and cj) : ^ J2 is a difjeomorphism, then 
the following diagram commute: 



method 



{{qk,Pk)}k 



method 



{{3k,rk)}k 



Proof. For simplic ity we d enote the map T*(f) : T*£2 T*y hy x- Further, we 
denote Algorithms ISJp^s] on T*^ for the triple (g, d,V)hy $^g, and $rks- 



The corresponding algorithms on T*S^ for the pulled back triple (h,e, H^) 
{(l>*g,(j)*d,4>*V) are denoted Tgg, Tjg and Tj^j^g. We need to show that 4>3g 



3S°X 



^2S 



2S°X 



-1 and $^j,s 



XoT 



RKS ° X" 



First, notice that the vector fields Xt, sY and Xy on T*J2 are defined 
respectively by the triples (g, 0,0), (0,d,0) and {0,0, V). Thus, from the com- 
mutative diagram Q it follows that exp{tXu) = xo exp{tXT) ox~^ , exTp{tY) = 
X o exp{tZ) o x^^ and exp(tXv) = X ° exp(tXw) o x~^7 where Xjj, eZ and 
X\Y are the vector fields on T*S' defined respectively by the triples (h,0,0), 
(0,8,0) and (0,0, W^). From the definition of Algorithm |3.1| it now follows 

that $(,'g ^xoT^s^X-'- 
Likewise, the vector field Xy + sY is defined by the triple (0, d, V). Thus, it 
follows again from diagram Q that cxpiUXy+eY)) = xoexp{t{Xw+£Z))ox^^ , 



3.2 



we get = X ° T^g o 



and from the definition of Algorithm 

Lastly, notice that Xy + eY has trivial dynamics in q. Thus, it reduces to a 
system on T* Q (only the p variable is affected by its flow) . From diagram ([t]) 
it follows that X\y + sZ corresponds to a change of variables q = (f>{s) and 
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r = T*(f> ■ p. This is a linear map in the nionientuni variable. It follows that 
"^Xv+eY — X° '^Xw+eZ ° since is a linear method, i.e., invariant under 
a linear change of variables (every Runge-Kutta method i s). I t now follows that 
*^RKS — X° "^RKS ° X^^ from the definition of Algorithm 3.3 □ 



Remark 3.2. Although invariance with respect to choice of configuration co- 
ordinates for numerical integration of system Q seems to be a highly natural 
requirement, most families of integrators do not have this property. For ex- 
ample, symplectic Runge-Kutta methods, which are invariant under a general 
linear change of coordinates and preserve the canonical symplectic form, are not 
invariant under co-tangent lifted point transformations. In engineering applica- 
tion, this has an impact, since the result of a numerical simulation then depends 
on whether Cartesian, cylindrical or spherical coordinates are used. 



4 Backward Error Analysis 

In this section we present a backward error analysis for the methods presented 
in Section |3] We start off with a brief review of the general notion of backward 
error analysis for ordinary differential equations. Thereafter, we show how the 
framework can be used for the class of Rayleigh damped problems considered 
in this paper. 

4.1 Review of Framework 

Let ^ he a, phase space manifold and X € X{^) a smooth vector field. The 
set of diffeomorphisms on ^ is denoted Diff(^) (this set forms an infinite 
dimensional group under compositions). Following previous authors |23 } I2 H [9]. 
we now define rigorously what is meant by an integrator. 

Definition 4.1. An integrator for a vector field X E X(^) is a one-parameter 
family : ^ ^ of diffeomorphisms that is smooth in h and satisfies: 
$° = Id and $''(z) - exp(/iX)(z) = 0{h''+^) for all z e ^, where r > 1 is the 
order of the integrator. 

Given an integrator for X £ X{^), the basic notion of backward error 
analysis is to find a modified vector field Xh, depending smoothly on h, such that 
its fiow coincides with i.e., such that exp{hXh) = Thus, the question of 
finding a modified vector field is equivalent to the question of finding an inverse 
of the exponential map exp : — > Diff(^^). However, it is well known that 
exp is not surjective, not even in an arbitrary small neighborhood of the identity 
(see e.g. [5]). Thus, consider instead the restriction to the real analytic case. 
Indeed, let ^ be a real analytic manifold (cf. [13]), let X°(^) C X(^) be the 
space of real analytic vector fields on and let Diff"(^) C Diff(^) be the 
corresponding real analytic diffeomorphisms. Even in this case, the restricted 
exponential map exp : — > DiS°'{^) is not locally surjective |13l Chap. IX]. 
However, one can still find a formal modified vector field by asymptotic expan- 
sion in the step size parameter: 

Xh = Xo + hXi + h^X2 + ... (9) 
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where Xo,Xi,X2,... G X"(^). Given an integrator for X e the 
vector fields Xo , , X2 , . . . are defined recursively by 

= iTo 

where Xq = X and X^^k ~ Y^^^a^^^i- Thus, with this construction <^^{z) — 
exp{hXfi^k){z) — 0(/i''+^), which follows from Taylor's theorem. Notice that if 
the integrator is of order r, then Xk — for < fc < r. 

In general the formal series ^ does not converge, so instead one has to 
truncate the series. The basic result in backward error analysis states that there 
is an optimal truncation index k, depending on h, such that — cxp{hXfi^k) is 
exponentially small. There are many contributors to this result: we refer to |5J 
Chap. IX] and references therein for an account of its origin. The original result 
is obtained for Euclidean phase space ^ = W . For the purpose of this paper, 
we use a generalization to manifolds given in [5] : 

Theorem 4.1. Let ^ be a real analytic manifold, \J G ^ a compact subset, 
X G X''(^^) a real analytic vector field, and an integrator for X such that 
h I— )• <i>''(z) is real analytic for all z ^ U. Then there exists a distance function 
dist : X M+, a truncation index k (depending on h), and constants 

C, 7, /iQ > such that 

dist{^''{z),exp{hXh,k){z)) < hCe-^/'' 

for all z whenever h < Hq. 



Remark 4.1. The distance function used in Theorem |4.1| is obtained by first 
embedding 13'' in (by the Whitney embedding theorem) and then restricting 
the Euclidean distance function on M" to the embedded submanifold. In par- 
ticular, in the case when ^ = M" , the distance function used is just the usual 
Euclidean distance. See [9] for details. 

Let us now consider the special case when ^ = T*^ (as in this paper). 
Denote by XHam(r*i?) the space of Hamiltonian vector fields on T*i2 (with 
respect to the canonical symplectic form), i.e., X £ XHam(r*i?) implies that 
X = Xh for some Hamiltonian function H £ J'{T*JS). Correspondingly, we 
have the subgroup of exact symplectic diffeomorphisms DiffHam(2^*^)- It is 
well known that exp(X) e DiffHam(r*^) if and only if X G XHam(r*^). 

Remark 4.2. Notice that there is, in general, a difference between Hamiltonian 
vector fields XHam(2^*=S) and symplectic vector fields Xsp{T*^), the former be- 
ing a subalgebra of the latter. Indeed, X G Xsp{T*^) means that X preserves 
the symplectic form, whereas X G XHam(T*^) means that there exists a glob- 
ally defined H G S'(T*^) such that X — Xh- Correspondingly, the group of 
exact symplectic maps DiffHam(T'*^) is a subgroup of the group of symplectic 
maps Diffsp(T*^). See e.g. [131122] for further reading. 

Now, let be an integrator for Xh G XHam(r*^). Naturally, is called 
a symplectic integrator if <i>'' G Diffgp (T*i2) for each fixed h > 0, and it is called 
an exact symplectic integrator if G DiffHam(r*i?) for each fixed h > 0. In ad- 
dition to Theorem |4.1[ another basic result in backward error analysis is that if 
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is an exact symplectic integrator, then each of the vector fields in the formal 
series ^ belong to XuamiT*^). Thus, each of these vector fields correspond 
to a Hamiltonian function. For a real analytic Hamiltonian system integrated 
by an exact symplectic integrator, application of Theorem |4 . 1 1 yields that if the 
numerical solution stays on a compact subset, then there exists a truncated mod- 
ified Hamiltonian function (corresponding to a truncated modified Hamiltonian 
vector field) which is conserved by the numerical solution up to an exponentially 
small term for exponentially long times (see e.g. [B] Sect. IX.8]). Furthermore, it 
can be shown that the family of truncated modified Hamiltonian functions have 
a common Lipschitz constant, independent of the step size |H1 Sect. IX. 7. 2]. 

A deeper consequence of the modified vector field preserving the Hamiltonian 
structure is that if Xh is Arnold- Liouville integrahle (cf. [HI Sect. X.l]), i.e., it 
has invariant tori, then an exact symplectic integrator will preserve perturbed 
tori for exponentially long times (under certain technical assumptions associated 
with KAM theory). 



4.2 Application to Perturbed Hamiltonian Problems 

In this section we derive the modified vector fields for the numerical methods 
suggested in Section [3j Our aim is to study the energy dissipation behavior of 
the methods, using a backward error analysis approach. The idea is to study the 
evolution in the numerical solution of a truncated modified Hamiltonian func- 
tion. In particular, we show that up to an exponentially small term, the sug- 
gested methods have asymptotically correct energy behavior in the dissipation 
parameter e. Furthermore, the evolution of a truncated modified Hamiltonian 
for Algorithm |3.1[ is monotone up to an exponentially small term. 

For convenience, we say that a numerical method for perturbed Hamiltonian 
systems of the form Q is (exact) symplectic, if it is (exact) symplectic for e = 0. 
Thus, the suggested methods Algorithm |3.1f|3.3| are exact symplectic, which 
follows directly from the splitting approach, and the fact that DiffHam('7^*=S) 
forms a group under composition. Throughout this section, we assume that 
phase space T*=S carries a real analytic structure. 

Let = Xt + Xv + sY be the vector field in equation [8] From equation [s] 
it follows that 

77(exp(tZ,)(z)) - H{z) = / d,(,)(p(s)«,p(s)«) ds (10) 



where (g(s),p(s)) = exp(sZe)(z). Our aim is to find a modified analog of (10) 
corresponding a numerical integrator, ideally by replacing the Hamiltonian with 
a modified Hamiltonian, and the dissipation tensor with a modified dissipation 



tensor. In particular, there are two qualitative features of equation (10 1 that we 
would like to find analogs of in the numerical solution: 

• iJ (exp(tZj)(z)) — H{z) is proportional to e; 

• H{exp{tZ^){z)) ~ H{z) < for i > 0. 

We first give the general structure of the formal modified vector field corre- 
sponding to a symplectic method: 
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Lemma 4.2. Let be an exact symplectic integrator of order r for problem Q . 
Then its formal modified vector field is of the form 

oo 
k=r 

where Xh^. G XHam(T*=S) and Yk^e G X(r*=S) depending smoothly on e. 
Proof. For £ = the formal series is of the form 

oo 

Zhfl = Xt + Xy + ^ H'^Xh^ , 
k—r 

since the method is exact symplectic. The result now follows from Taylor's 
theorem since depends smoothly on e. □ 



From this result we get a modified analog of equation (10 1 for general exact 
symplectic integrators. Indeed, let H^.n = T + V + J2k=r ^'^^k, with Hk as in 
Lemma |4.2[ Then we have the following result: 



Theorem 4.3. Let be an exact symplectic integrator of order r for prob- 
lem and let U C T*^ be a compact subset. Assume that = Xx+Xy+sY 
is real analytic in U, and that h i— > is real analytic for every z G U. Then 

there is a truncation index N depending on h and e such that 

HhA'^'^iz)) - Hu,n{^) + e / dh,N,e{eMsZh,N,e){z)) ds < /lACe-^/'' 

whenever h < Hq and e < Sq, where A, C, 7, /iq, eo > are constants not depend- 
ing on h or e, and where 

N 
k=r 



for Yk^e i'n Lemma 4-2 



Proof. From truncation of the modified vector field in Lemma [4?2| it follows that 

= {dHh,N{z),XH^^^(z) +£YhMAz)) 
= £{dHh,N(z),Yh^N.A^)) 

N 

= -e(d,(p«,p«) - ^ h\dHh,N{z). n-,s(2)> ). 

fe=r 



Thus, 

ft 



Hh,N{exp{tZh^N,e)iz)) - Hh,N{z) = / dh^N,e{(i'y^V{sZh,N,e){z)) d 



is. 
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Moreover, we have 

Hh,N{'^>''{z)) - Hh,Niz) = Hh,N{exp{tZh,N.,e){z)) - Hh,N{z) 

+ Hh,N{^''{z)) ~ Hh,N{<iMtZh,N,e){z)), 

so we get 

\Hh,N{^^{z)) ~ Hh,N{z) + e I 6h,N.e{Qy^V'{sZh,N.e){z))<^s\ 

Jo 

= \Hh,N{^''{z)) - HhM{<^MtZh,N^e){z))\ 

Let A > be an /i-independent Lipschitz constant for H^^n with respect to 



the distance function dist in Theorem 4.1 (such a Lipschitz bound is known to 



exists for truncated modified Hamiltonians, see e.g. [8, Theorem 8.1, Sect. IX. 8]). 
Thus, an estimate for the right hand side is 

\Hh,N{^''{z)) - Hh.N{eMnhZh,N.e){z))\ < A dist($''(2), exp(/tZ,,,jv.e)(2)) • 

(11) 

Now, in order to get a bound which is uniform in e, let eo > and consider the 
extended map e) := e) on U x [0, Eq] (notice that depends on e), 

which is an integrator for the extended vector field Z{z,e) := (Z^{z),Q). Since 
the extended e-part is integrated exactly by it holds that exp{hZii j^){z, e) = 
(exjp{hZh^N.e){z),s), so we have 

dist($''(z,e),exp(/ii-,i,7v)(2,£)) = (dist($''(z), exp(;iZft,Ar^e)(z))^ + |e - e]' 



Thus, in combination with (11 1 we have 



\Hh,N{^''{z)) - Hh,N{exp{nhZh,N,6){z))\ < X^t{^''{z,e),exp{hZh,N){z,e)) 



for h < ho > 0, where the last estimate follows from Theorem |4.1| applied to 
the extended integrator which can be done since U x [0,eo] is compact and 
Z{z,e) is real analytic on U x [0,eo]- This proves the theorem. □ 

From the result we see that, up to an exponentially small term, the modi- 
fied energy evolves with a rate which is 0{eh^) close to the energy dissipation 
rate of the exact solution. Thus, the dynamics is asymptotically correct in the 
perturbation parameter e. 

In transient regions of the phase space, i.e., regions where the value of dissi- 
pation tensor d is relatively large. Theorem |4 . 3 1 asserts that the energy evolves in 
a monotone fashion. More precisely, let U C r*=S be a compact subset such that 
sup^^u\(ih.N,e{z) — dq{p^,p^)\/dq{p^,p^) <^ 1 for the chosen step size. As long as 
the numerical solution stays in U, Theorem |4.3| then asserts that the modified 
Hamiltonian is decreasing from step to step, and that the numerical dissipation 
rate is accurate relative to the exact dissipation rate. Thus, the transient phase 
in the simulation is captured very well with a symplectic integrator, which is 
often important in application, for example when accurate estimates of energy 
losses are essential. 
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After the transient phase some type of "steady state" dynamics is commonly 
reached. In phase space, this often occur close to the subset kerd = {z E 
T*^;dq{p\p^) = 0}, i.e., where the dissipation vanishes. Furthermore, it is 
common that dissipative forces (e.g. friction forces) are small relative to conser- 
vative forces. In such regions, one does not necessarily have that Ahjs[^g{z) > 0, 
so Theorem 14.31 does not assert that the evolution of the modified Hamiltonian 
is monotone. However, for the special case of the suggested Algorithm |3.1[ we 
now derive a refined version of Theorem |4.3[ for which monotone decrease of 
the modified Hamiltonian is asserted for small enough step sizes (up to an ex- 
ponentially small correction term). The result is based on the following two 
lemmas: 

Lemma 4.4. The formal modified Hamiltonian for the Stormer- Verlet method 
'^^y applied to problem (jsj with e = is of the form 



OO oo 



Y^Y.h'''^'-'^{Z2,.2,UpK...J) 
k=i 1=1 ^T! ' 

Ik copies 

where Vie € 3^{T*J2) and gik.ie € '^ki'^) '^'"^ totally symmetric tensor fields. 

Proof. If a € '^(■S) and b e 7^(=S) are two totally symmetric tensor fields, then 
the Poisson bracket between the functions {q,p) i— >■ aq{p^ , . . . , p^) and {q,p) i— 
bq{p^ , . . . , p"^) is given by (g,_p) Cq{p^ , . . . ,p^), where c S T2_|_£(^) is the totally 
symmetric tensor 

Cq{-, ...,•)= aq{dqhq{-, ...,■),-,...,■) + ...+ a,(-, dqhq{-, ...,•)) 

~hq{dqaq{-, ...,•),•,•■•,•) + ••• + bg(-, dqaq{-, •)). 

Since Stormer- Verlet is a splitting method, and since V and g are tensors, it 
follows from the Baker-Campbell-Hausdorff formula that the modified Hamil- 
tonian is a sum of terms of the form /i*^-'gij(p'*, ■ . ■ ,P^), where gij £ 7^{^). 
Further, since the Stormer- Verlet method preserves reversibility it holds that 
Hh{q,p) — Hh{q,—p). Thus, only the even order tensors remain. Also, only 
the even powers of h survive, since the method is symmetric. This yields the 
result. □ 

Lemma 4.5. Let Y £ X(T*=S) be the dissipative vector field in equation (Kj), and 



let Hfi N be the truncation of the formal modified Hamiltonian in Lemma 4-2 



Further, let U C T*^ be a compact subset. Then there exists Hq > such that 

£YHhM = -dqiEh,Niz)pKp^) < 

for all z £ U whenever h < Hq, where E}i^n{z) '■ Tq^ — >■ T^cS is an opera- 
tor which is self-adjoint and positive definite with respect to d, and such that 
EhM{z) - Id = 0(/i2) for all z £D. 

Proof Let ((•,•)) denote the metric g,(-,-)- Let C{q) : Tq^ Tq^ be the 
positive semi-definite self-adjoint operator defined by {{C{q)u,v)) — dq(u,v). In 
other words, C{q)u = dq{u, •)". 
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It holds that 

£Y{g2k,2l)q{p\ ■ ■ ■ = {g2k.2t),{C{q)pKp\ . . . + 

{g2k,2d,{pKC{q)pKp\...J)+...+ 
(g2fe,2£),(p», . . . ,p\C{q)p^) {A2k.2i),{pK ■ . ■ 

which defines the totally symmetric tensor d2fc,2£ G 'J^^(T*J2). Notice that the 
kernel of d2fe,2« is contained in the kernel of C{q), i.e., if it S ker(C(g)), then 
{<i2k,2e)qiu, ■) = 0. Next, we define the self-adjoint operator C2k,2eiz) : 
Tg^ Tq^ by ((C2fe,2^(2)u,w)) = {d2k,2i)q{pK ■ ■ ■ ,pKu,v). Again, notice that 
ker(C2fc,2^ (-z)) C ker(C((7)). All together, summing up the terms up to order 
in Lemma |4.4[ we get 

£YHnM{q,p) = -{{C{q)p\p^)) - {{Ct,r^{z)pKp^)) 

where C^j^{z) is a self-adjoint operator Tq^ Tq^, with ker(C^^(z)) C 
ker{C{q)). Since both C{q) and C^j^{z) are self-adjoint it holds that (ker C(q))^ 
is an invariant subspace for both of them (follows from the spectral theorem). 
Since d has constant rank and since U is compact, it follows that 

inf ||C(9)|(kcrC(g))^lP > 0. 

Thus, there is an Tiq > such that C{q) + h'^C^j^{z) stays positive semi-definite 
when h < ho and z S U; the zero eigenvalues stays zero, and the perturbation, 
acting in (ker C{q))-^, is small enough for its positive eigenvalues to stay positive. 
Hence, 

£YHh,N{z) = -{{{Ciq) + h'Ct,N{z))p\p^ < 

if h < ho and z e U. Next, define Eh,N{z) ■ Tq£2 TqJ2 by dq{Ek.^N{z)u,v) = 
{{{C{q) + h^Ctj^{z))u,v)) for u,v e (ker C{q))^ and Eh^N{z)kerC{q) = M. □ 

We now give the refined result on evolution of the modified Hamiltonian 
function for Algorithm |3.1| 



Theorem 4.6. Let U C T*J^ be a compact subset, and assume that g, d, V 



3.1 



there 
'dex N 



real analytic on U. Then, for the integrator defined by Algorithm 
exists ho > 0, depending on U,g,d,V^ but not on e, and a truncation in 
depending on h, such that 

rh/2 

,Ar(z)+£ / d/i.Ar(exp(s£F)(z)) ds 
Jo 

h/2 

d,,.w(cxp(s£y) o $^^^o exp(/i£y/2)(z)) ds\ < hXCe-^'^' 

whenever h < ho and z G U, where A,C, 7,/io > are constants not depending 
on £, and where 

dh^Niz) :=dqiEh,Niz)pKp^) > 



with Eh^j^{z) : TqJS Tq£2 as in Lemma 4-5 
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Phase diagrams for damped pendulum using Algorithms |3.1|43.3| 



h = 0.3 




Figure 1: Phase diagrams for the damped pendulum problems integrated with 
the suggested methods (they all overlap at this level of detail) for various step 
sizes. Notice that the amount of dissipation does not depend much on the step 
size (the size of the hole in the middle is roughly the same). 

Proof. We have 

Hh,N{^3siz)) ~ Hh,N{z) = Hh,N{exp{heY/2)iz)) - Hh^iz) 

+Hh,N i'^SY ° exp(tey/2)(z)) ~ Hh,N{exp{heY/2){z)) 

+HH,N{eMheY/2) o o cxp(tey/2)(z)) - i?h.,w($gv ° exp(ter/2)(z)) 

From Theorem 14.31 with e = it follows that 

Hh^i'^SY ° exp(ter/2)(z)) - HhM{eMheYl2){z)) < hXCe-^/'' 

whenever h < h'^ for A,C, 7,ft.Q independent of e (since it is zero when we 
apply the theorem). Since exp{heY/2) is the exact ft,-flow of eY/2 it follows 
from Lemma |4.5| that there exists Hq > (clearly independent of e) such that 
<ih,N{z) = dq{Eii^Niz)p^ ,p^) > 0. This gives the four remaining terms in the 
sum. Now chose ho = min(/iQ, Hq ). □ 

Notice in particular that this result asserts, up to an exponentially small rest 
term, that the modified Hamiltonian decreases monotonely for small enough step 
sizes. The result is optimal, in the sense that the same exponential term would 
remain if we take e = 0. 

4.3 Numerical Example: 2— D Pendulum 

Let ^ — R, gq{u,v) — ^uv, d = g and V{q) = 1 — cos{q). The system so 
obtained describes a damped non-linear pendulum in the vertical plane, with 



16 



Phase diagrams for damped pendulum using Heun's method 

h = 0.3 




Figure 2: Phase diagrams for the damped pendulum problems integrated with 
Heun's method for various step sizes. Notice that the amount of dissipation 
(corresponding to the size of the hole in the middle) strongly depends on the 
step size. 



17 



Evolution of energy error {h = 0.2) 




Comparison of energy in between methods {h = 0.2) 

10-^ 5 . . . . 




(Upper) Evolution of energy for Algorithms 3.1 3.3 



and for Heun's 
whereas it grows 



Figure 3: 

method. The error is bounded for the suggested methods, 
linearly for Heun's method. (Lower) Comparison of energy evolution in between 
Algorithms 3.1 3.3 The difference is about a factor lO'^ or 10'* smaller than the 



actual energy error. 
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unit mass and unit gravity (see e.g. [5J Sect. I.l] for details on this problem, in 
the un-damped case). 



The system is integrated with each of Algorithm 3.1 3.3 as well as with 



Heun's explicit second order method. (We also use Heun's method as the choice 



of Runge-Kutta method for Algorithm 3.3 ) The following data was used: q(0) 



0.9 TT, p(0) = 0, £ = 10~ and integration time interval [0,200], for various step 
sizes h e [10~\ 1]. 

Phase diagrams for various step sizes using Algorithms |3.1f|3.3| are given in 
Figure [l] There is virtually no difference between the suggested methods, i.e., 
the phase curves overlap. Notice that the size of the "hole" in the middle of the 
phase diagrams is almost independent of the step size. This verifies the result in 
Theorem |4.3| that the dissipation rate is asymptotically correct independent of 
the step size. The corresponding phase diagrams for Heun's method (sometimes 
also called Runge's method|^are shown in Figure [2j Notice that the dissipation 
rate in all the simulations is to small, and depends heavily on the step size. This 
reflects the fact that Heun's method (as most explicit Runge-Kutta methods) 
assembles "numerical" energy. 

A detailed plot of the energy error as a function of time, is shown in the 
upper graph of Figure |3] (The "correct" energy is computed by highly accu- 
rate numerical integration.) In particular, notice that the energy error for the 
suggested methods seems to stay bounded, whereas it grows linearly for Heun's 
method. Although this result is not fully explained by our backward error anal- 
ysis, it is related to the result obtained in Theorem |4.3| and Theorem |4.6| that 
the dissipation rate for the modified Hamiltonian is asymptotically correct in e. 

The lower graph in Figure [3] contains comparison of the energy between the 
suggested methods, i.e., 

|i/o<i>^g-i/o$^j,sl' |i?°<i>Ss-^°*2sl' I^^°*3S-^°<KSI 

(this difference is too detailed to recognize in the upper graph) . Comparing with 
the values in the upper graph, notice that the difference in computed energy 
between the methods is much smaller than the actual energy error, i.e., they 
produce practically the same result. Also, notice that among the methods, the 
difference in result between $23 ^^'^ ^rks indistinguishable as compared to 
<i>3g, i.e., <i>3g is the "out-her" of the three. 



5 Symmetry and Conservation of Momentum 

In this section we show that the methods suggested in Section |3] preserve in- 
variance under a symmetry, and conserve corresponding momentum maps. For 
preliminaries on geometric mechanics, see |17] or [T]. 

Let be a Lie group acting on ^ by : x ^ — > We write <f)g := (j){g, •). 
Further, for 5 S we use the notation g ■ q = 4>g{q) and g ■ z — <j)^,g{z) :— 
{T*,^(j)g-i){z) for the lifted action. The action is associated with a corresponding 
momentum m,ap J : T*£2 — q* . where g* is the dual of the Lie algebra g of . 
Explicitly, the momentum map is defined by 

(J(q,p),e) = (p,ei2(9)>, veeg, (12) 








^ The Butcher tableau is 1 


1 




1/2 1/2 
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where is the infinitesimal generator of the action on i.e., £,s{q) '■— 

^|t=o(exp(tO-g). 

Now, let Y G X(r*=S) and assume that its flow is -equivariant, i.e., that it 
commutes with the action: 

exp(iy) o (f)^,g — (j)^,g o exp(ty). (13) 

The infinitesimal version of this relation is [i^, ^t*^] — for all C G fli where 
Ct*^(-z) '■— -37|t=o(exp(i^)-z) is the infinitesimal generator of the action on 
The vector field Y is then called '^^ -invariant. Since the action on r*=S is 
canonical, the vector field ^t*^ is symplectic. In fact, it corresponds to the 
Hamiltonian function J(^) = (J(')' 0' i-^-i Ct*^ = ^J(4) so it is exact symplectic 
(see [If, Sect. 11.2]). 

Associated with the momentum map is: 

• The subgroup of $^-equivariant diffeomorphisms 

Diff^^(T*^) = {ipe Difr(T*^); o = o 95, Vg e 

Since ^t*^ generates 0*oxp(t5)j the corresponding subalgebra of vector 
fields is given by 



• The subgroup of momentum conserving diffeomorphisms 

Diffj(T*^) = e Diff(r*^); j(0 o ^ = e 0}, 

with corresponding subalgebra of vector fields given by 

Xj(r*^) = {xg X(r*i2); £xJiO = 0, e g}. 

Notice that ^^-equivariance does not necessarily imply conservation of the 
momentum map, not even in the Hamiltonian case|^ Nor does conservation 
of momentum imply ^^-equi variance. However, if a function H G 3'(T*^) is 
^^-invariant, i.e., {H, J(^)} = for all S, E Q, then for the associated Hamilto- 
nian vector field it holds that Xh € X<^iT*^) n Xj(T*^) since [Xn^^T'd = 

X{H,Jim = and £x« J(0 = {H, JiO} = 0. 

Let X<^,j(T*^) := X^^(r*^)nXj(r*^). We now give criterions for the vector 
fields involved in system ([s]) to be ?^-invariant and/or momentum preserving. 

Proposition 5.1. Let g, d and V be the tensors used in equation and Xt, 
Y and Xy the corresponding vector fields in equation (|8|. Then, 

£i^g = =^ XTeX<^j{T*^) 
£i^V = =^ XveX^jiT*^) 
i^^d = =^ y e Xj(r*^) 



^For example, take S = = R, H(q,p) = q and J{Oi<l,P) = P ■ 5- Then [Xh,Ct*.s] 
= = 0, but £x„J{0 = {H,J{0} = C 7^ 0. 
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Proof. The result is already well known in the Hamiltonian case of Xt and 
Xv- Indeed, £^c^g — implies that acts on T*^ by isometrics. In turn, this 
implies that T is ^^-invariant. Further, i^j^l^ = implies that V is ^^-invariant. 
In turn, Xt and Xy are ?f-invariant and {T, J($)} = {V, J{^)} = 0. (See |24] 
or P Sect. 4.5] or jTH Chap. 3] for details). 

The result remains to be shown for the non-Hamiltonian vector field Y . First 
the third implication. Since Y{q,p) — (0, — dg(p', •)) and J{£,){q,p) = {p,S,.s{q)), 
it holds that £yj(0(g,p) = -dgipK^Jsiq))- Now, d,(p», ^^(g)) = (i5^d),(p«) = 
by the presumption. Thus, £yJ(C) = 0, so y G Xj{T*^). 

Next, the last implication. Let f denote the bi-hnear form T*^ x T*=S — > M 
given by f,(p,r) := g,(p«,r«) = (p,r»). If i^^^g = then = (see \m 

Sect. 3.1]). A calculation equivalent to the one in equation ^ yields 

(^,,),y(g,p) = (o,-(<^,gd),((<^,/),(p, •),•))• 

Chose now the path g{t) — t^. replace g for g{t). and take the time derivative 
-g^lt^g. The left hand side then becomes the Lie derivative of Y in the direction 
^7^*^. Using the product rule on the right hand side we get 

£^^.^Y{q,p) = (O, -(£e.d),(f,(p, •), •) - d,{{£^J),{p, •), •))• 

By the presumption the right hand side vanishes, so [^T'^,y] = 0. Thus, 
Y e Xgp(T*^). □ 

We now have the following basic result, which concerns conservation of mo- 
mentum and ^^-invariance for the suggested methods. 

Proposition 5.2. Let denote any ofBiSr^iT*^), Diffj(T*^) or Diff^^jT*^, 
and let [) denote the corresponding suhalgebra of vector fiel ds. Then Xt,Y,Xv G 
() implies that ^^g, and ^'jij^g, defined by Algorithms 3.T 3.3 belong to . 



Proof. The result is clear for Algorithm |3.1}|3T2| since is closed under com- 
position. Moreover, for $p^j^g = exp(/iX7'/2) o "^^^^^y o exp(/iXr/2), it holds 
that the exp(/iArT/2) € if Xt € \). The vector field Y defines a linear 
system on T*^. Since the momentum map is linear in p, and since all Runge- 
Kutta methods conserve linear invariants, it holds that ^'^^_|_gy € Diffj(r*.S) if 
Xv,Y e Xi{T*^). Furthermore, we notice that t he lif ted action is exactly 
the co-tangent lift of 4>g-i- Thus, by Proposition 



3.1 



*^RKS invariant under 



the change of coordinates {q,p) = 4i^g{s,r). That is, 

eMH4>*a)*XT/2) o ° eMH^*g)*XT/'2) o c/j.g = o ^j^^S- 

Now, if Xt.Y,Xv e Xc^iT*^), then (0*g)*XT = Xt, {<l)*g)*Y = Y and 
(0*g)*Xy = Xv, so $^j^g o (j)^g = o $j^j^g, i.e., $^j^g e Diff^(r*^). □ 



In particular, if g, d and V fulfill the requirements of Proposition |5.H then 
it follows by Proposition |5.2| that Algorithms |3.1f|3.3| are both ^-equivariant 
and momentum conserving. It is important to point out that this result does 
not hold for general exact symplectic integrators, not even if the integrator 
is momentum preserving when s = 0. For example the alternative splitting 
method 

*RKS-B = eMhXv/2) o vl/^^^^^ o eMhXv/2) (14) 
reduces to Stormer-Verlet when e = 0. but does not conserve momentum. 
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6 Near Conservation of Relative Equilibria 



In this section we investigate how the methods suggested in Section [3] preserve 
relative equihbria. For background on geometric reduction theory see e.g. |15l 
Chap. 3], Pa Chap. 4], or [H Chap. 1-2]. 

Let ^ be a Lie group acting on =S as in the previous section. We assume 
that J# acts freely and properly on (cf. [151 Chap. 3]), which roughly speaking 
means that the action is non-singular. 

Definition 6.1. Let X e X{T*J2). A solution curve j{t) to z — X{z) is called 
a relative equilibrium if there exists a ^ G such that ^{t) = exp(t^) • 7(0). 

The study of relative equilibria of Hamiltonian systems on a co-tangent bun- 
dle T*J2 is closely related to the theory of co-tangent bundle reduction. The 
theme is to "quotient away" the part of the dynamics generated by the sym- 
metry group and thus, for each /i € 0*, obtain a reduced phase space given 
by — J^^(/i)/G'^, where is the subgroup of that leaves fj, invariant 
under the co-adjoint action, i.e., = {g G Ad*-i/i — fi}. 

Remark 6.1. If is an Abelian Lie group (all elements commute), it holds 
that % = ^, so = J-i(^)/^, and dim = 2dim^ - 2dim^. In this 
case, the reduced phase space is isomorphic to T*(=S/?#), equipped with 
a non-canonical symplectic structure cj^, depending on the momentum map 
value ^. It holds that luq is the canonical symplectic structure on T* /'^). In 
general, = a;o + /3p, for a 2-form /3p usually called magnetic term. In classical 
mechanical examples, where typically is a rotation group, corresponds 
to centrifugal forces due to the rotation. If dim^^ = dim=S we have n first 
integrals in involutions, i.e., the system is Arnold-Liouville integrable. The 
modern notion of co-tangent reduction in the Abelian case was developed by 
Smale [24j . However, its roots goes back to Lagrange, Poisson, Jacobi, and 
RouthQ 

Let : J^^(/i) r*i? be the natural inclusion, and tt^ : J^^(/i) — s- 
the projection 7rp(z) = [z], where [z] G = J~^(/i)/^^^ is the equivalence 
class of z G J^^(/x). Given X G X^j(T*i?) and initial data on J~^(/i), the flow 
exp(iX) restricts to the invariant momentum manifold J^^(/i). Further, due to 
the symmetry, it drops to a flow exp(tAr) on Diff(,^3^^) that fulfills 7r^oexp(iX) = 
exp(iX) o TT^, where X is the vector field on defined by Ttt^ o X — X o tt^. 
Notice that for these constructions to make sense, it is essential that X is both 
^^-invariant and momentum conserving. If ^(t) is a relative equilibrium, then 
TT/i o 7(t) := [7(t)] = [7(0)]. Indeed, "f{t) is a relative equilibrium solution of the 
dynamical system i = X{z) if and only if [7(0)] is an equilibrium of [z] = 
i.e.,X([7(0)])=0. 

Note that if ^{t) = {q{t),p(t)) is a relative equilibrium for system ([s]), then 
it must hold that dq(()(p(i)'*,p(i)'*) ~ for all t, since the motion 7(t) conserves 
the Hamiltonian. 

Now, fix some e > and let Z = Xt + Xy + be the vector field for the 
system (|8]), and assume that g, d, V fulfill all the requirements in Proposition 5.1 



so that Z G X<^.j(r*=S). Further, assume that system ([8| has an asymptotically 



^ A historical account of the theory of co-tangent bundle reduction is presented in the 
introduction of the monograph by Marsden et. al. |16j . 
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stable relative equilibrium j{t), i.e., the critical point [7(0)] of Z is asymptoti- 
cally stable. In general, an integrator for Z does not have a nearby relative 
equilibrium, not even for very small step sizes. The reason for this is that Z is 
not structurally stable with respect to arbitrary perturbations in X(T*JS). That 
is, for a small perturbation eZ^ £ X{T*^) the vector field Z + eZ^ does not, 
in general, possess a perturbed relative equilibrium. However, the stationary 
point [7(0)] of the reduced vector field Z is structurally stable with respect to 

perturbations in X{^fj,). That is, if eZ E X{3^^), then for small enough e, the 

perturbed reduced vector field Z-\-eZ^ has a nearby stationary point. As a con- 
sequence, the relative equilibrium 7(t) of Z is structurally stable with respect 
to perturbations eZ^ £ Xcg j{T*^), since then the vector field Z + eZ^ drops to 

Z + eZ^ under symmetry reduction. The corresponding argument also works 
on the level of flow maps. Indeed, let e Diffc^ j(T*^) be a ^I'-equivariant 
and momentum preserving integrator, approximating the exact flow exp(/iZ). 

By the quotient map, the integrator drops to a map $ G Diff(^^), which 
approximates exp(/i.Z) S Diff(^^), and since the fix point [7(0)] of exp(/iZ) 

is hyperbolic, the perturbed map has a nearby fix point [7^(0)] for h small 
enough. In summary, we have the following: 

Proposition 6.1. Let g, d, fulfill all the requirements in Proposition \5.1\ and 
let he a -equivariant and momentum preserving integrator for system^ 
Assume that j{t) is an asymptotically stable relative equilibrium. Then pre- 
serves a nearby relative equilibrium for small enough step sizes. 

Notice that the above result heavily relies on hyperbolicity, which in turn 
requires that the dissipation parameter e is strictly positive and the reduced 
dissipation tensor d is non-degenerate, since an equilibrium of a Hamiltonian 
system (typically) is elliptic. Thus, the requirement is that the step size is small 
enough in relation to the perturbation eY, i.e., essentially that ^ e. 

Based on Theorem |4.6| we now derive a refined result, which states that 
Algorithm |3.1| preserves a stable nearby relative equilibrium for exponentially 
long times, independent of e. Recall that a relative equilibrium ^(t) is stable if it 
corresponds to some local minimum [7(0)] of the reduced Hamiltonian H (that 
is, dH{[j{0)]) = and the Hessian of H at [7(0)] is strictly positive definite). 

Theorem 6.2. Let g,d,y be real analytic and fulfill all the requirements in 



Proposition 5.1 Assume that ^{t) is a relative equilibrium of system ([8]), such 
that the critical point [7(0)] of Z is stable, and let U be any compact neighborhood 
of [7(0)] such that the reduced Hamiltonian H is strictly convex on U. Then 
there exists another neighborhood \/ d \J of [7(0)] and constants Hq^k^K > 
0, independent of e, s uch t hat the reduced numerical solution [zq], [zi], . . . ,[z„ 



generated by Algorithm 3.1 stays in U over exponentially long times nh < Ke'^/^ 



whenever [zq] G V, h < Hq and e > 0. 

Proof. Since the integrator <i>3g defined by Algorithm 3.1 is ^^-equivariant and 



momentum preserving, it drops to a reduced integrator $33'' £ Diff(^p). Fur- 
ther, there is a reduced modified Hamiltonian given by Hh,N{[z]) ■— Hh^n{z) 
for any z £ [z] (well defined since iJ^.w is ?#-invariant). Since H is strictly 
convex on U, we may chose /iq > such that Hh,N is also strictly convex on 
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U whenever h < h'^, thus attaining a minimum Cmin = Hh.N{[z*i^ jvl) ™ 
terior of U. Next, let Cmax = inf [^jgog ^^/i.w ([2]). It holds that c^ax > Cmin- 
Now, let Cijiid = (cmax " Cmi„)/2. It holds that the set V = {[z] S U : 
Hh,N{[z\) — Cmid} is a neighborhood of [2,* jy] contained in the interior of U. Let 



[zo] G V. By Theorem 4.6 there exists /iq , 7, C, A > independent of e such that 
i^^^,^'([4'3s(•^)]) - Hh,N{[z]) < hCe-^/^ whenever [z] € U and /i < min(/i^, h'^)., 
which implies that i?;i.iv([z„]) — ^f;i.7v([^o]) ^ nhCe^'^/^ . Since [zq] S V it holds 
that i?/i,Ar([2n]) < Cmax whenever nhCe^'^^^ < Cmax ~ Cmid, which yields the 
result since Hh^N{[zn\) < c^ax implies that [2„] e U. □ 

The result states that if 7(t) is a stable relative equilibrium solution to sys- 
tem [8] then, for small enough step sizes, the numerical solution generated by 



Algorithm 3.1 will "almost" (for exponentially long times) preserve a modified 
relative equilibrium which is close to the exact one. In contrast to Proposi- 
tion [O] where the analysis is based on structural stability of hyperbolic critical 
points, the step size restriction in Theorem |6.2| is independent of e. In par- 
ticular, the result holds in the conservative case when e = 0. Furthermore, 
no non-degeneracy assumption on the dissipation tensor d is made: it may be 
degenerate in any direction. 

6.1 Relation to Symplectic Integration of Problems with 
Attracting Invariant Tori 

As mentioned in the introduction, there exists already a well developed theory 
for symplectic integration of Arnold-Liouville integrable systems, perturbed in 
such a way that only one invariant torus persist, and becomes attractive. For 
a thorough treatment of this theory, see |H1 Chap. XII]. In this section we give 
a short review of that framework, and compare it to the setting and analy- 
sis presented in the current paper. The two approaches are related, but not 
overlapping. 

Consider a perturbation of an integrable Hamiltonian system, which in 
action-angle variables (cf. [8, Chap. X]) can be written 

a = er(a, 6) 

aeM",6>eT". (15) 

= u>{a) +ep{a,e) 

Further, assume there is a point a* such that the frequencies uj{a*) are dio- 
phantine with exponent v (cf. [8, Chap. X]), and such that the angular aver- 
age r{a*) := r{a*,6)d0 is small and its Jacobian A = f'(a*) has strictly 
negative real part. Then system [15] has an invariant torus which attracts a 
neighborhood of {a*} x T" with exponential rate proportional to e. Now, [5J 
Theorem 5.2, Chap. XII] states that if a symplectic integrator is applied to 



system (Il5| expressed in canonical variables {q,p) E M^", then the numerical 
solution has a modified attractive invariant torus as long as /i'' < co|loge|~'^, 
where k = max(i^ -\- n -\- l^r) and Cq is a constant independent of h, e. 

Although Arnold-Liouville integrable systems are related to relative equilib- 



ria, in the sense that the former is a special case of the latter (see Remark 6.1 1 



systems of the form (15 1 are only overlapping with systems of the form i\8[ in 



the trivial case when the attractive manifold consists of only one single point 



24 



(corresponding to an asymptotically stable equilibrium). Indeed, the setting in 
this paper allows for a degenerate dissipation tensor, but has the requirement 
that the perturbation vector field eY fulfills £yH < for the Haniiltonian H. 
In contrast, the setting in j8j Chap. XII] does not require that the perturba- 
tion is monotonely decreasing energy, but instead that f in equation (15 1 is 
non-degenerate. As a consequence, in order for system (|8| to overlap with sys- 
tem (15) it is needed that the perturbed vector field fulfills fyH < 0, i.e., that 
the dissipation tensor d has full rank. In the common case of systems that have a 
locally minimal energy point in phase space, this means that the corresponding 
attractive invariant manifold is must be an asymptotically stable equilibrium 
point. 

As an example, take the Van der Pol equation studied in [8, Chap. XII]. 
Here T*^ = K^, with the harmonic oscillator Hamiltonian H{q,p) — p'^ /2+q'^ /2, 
and the perturbation is described by the tensor dq{u,v) = (1 — q^)uv (i.e., the 
perturbation is linear in p just as in this paper). However, for this tensor it 
does not hold that £yiJ < 0, since the tensor switches signature at q'^ — 
1. This is indeed the very reason that the Van der Pol system has a non- 
trivial attractive invariant manifold: the requirement £yH < would force 
the invariant manifold to be the equilibrium point {q,p) = (0,0). Likewise, the 
example considered in Section 6.2 below does not fit the setting in [F, Chap. XII], 
since the unperturbed system is not integrable, and since the dissipation tensor 
is degenerate. 

In conclusion, one can say that the main differences between the type of 
systems analyzed in |H1 Chap. XII] versus the ones in this paper, is the form of 
the perturbations, and that the former analysis requires the unperturbed system 
to be near integrable, whereas the latter requires the unperturbed system to be 
of standard form. 

Furthermore, the results of the actual numerical analysis is of different char- 
acter. This paper gives results on asymptotically correct dissipation rate (The- 
orem 4.3 and Theorem 4.6 1, and on near preservation of stable relative equilib- 



rium for exponentially long times independent of the perturbation parameter e 
(Theorem 6.2 1. The analysis in Chap. XII] gives results on near preservation 
of attractive invariant tori for indefinite time, but with step sizes depending 
weakly on e > 0. Also, as pointed out in [5J Chap. XII], any numerical integra- 
tor for system ( 15 ) preserve a nearby attractive manifold as long as h"^ ^ e. In 
contrast, for systems with relative equilibria the condition ^ e is not enough 
in order for a general integrator to preserve nearby relative equilibria. Indeed, 
as derived in Proposition 6.1 above, the integrator should be ^^-invariant and 
momentum preserving to assert that the numerical solution has nearby relative 
equilibria when <^ e. 



6.2 Numerical Example: Elastic 3— D Pendulum 

In this example we consider an elastic pendulum which is affected by gravity. 
The configuration space is ^ = M'^\{0}. There is a small Rayleigh damping 
term for the spring, but otherwise the pendulum is not damped. Thus, the 
dissipation is only active in the direction of the spring. The e parameter takes 
the role of the damping coefficient for the spring. Using Cartesian coordinates 
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the equations of motion are: 



= 

m 

where m is the mass, g e M'^ is the gravity, k is the stiffness of the spring, and 
£ is the length of the "unstretched" pendulum. Notice that the damping term 
only acts in the opposite direction of q. Also, when p is perpendicular to q the 
damping vanishes. In terms of the presented framework we have 



g,K.) = '^u'v, d,{u,v) = V{q) = \{i-\q\f q^g. 



Throughout the rest of the example, we work with unit constants: m — k — £ = 
1, and g = (0,0,-1). 

Remark 6.2. For e = the elastic pendulum problem is a non-integrable 
Hamiltonian system, meaning it has chaotic behavior; see Chap. 5]. 

The problem has an 5^-symmetry by rotation about the z-axis. Indeed, 
the setting is '^^ ~ with action on =S generated by ^^{q) — £, x q, where 
63 — (0,0,1). It is straightforward to check that g, d and V fulfill all the 



conditions in Proposition 5.1 Thus, the flow evolves on Diffc^j(T*^). The 



corresponding conserved momentum map is the azimutal angular momentum 



given by Jy,{q,p) = (q x p)3. In addition, it follows from Proposition 5.2 that 
Algorithms |3.1f|3.3| applied to this problem exactly conserves momentum. 

The problem is integrated with Algorithms |3. IffS . 3| and Heun's method, with 
the following initial data: q^ = (0, 1.55884573, -0.6) and = (1.34164079,0,0). 
The energy behavior for various step sizes and choice of e is shown in Figure |4] 
Notice that the qualitative behavior of Algorithms |3.1H3.3[ in terms of mean 
energy dissipation rate, is superior to the result with Heun's method. Indeed, 
contrary to Heun's method, the mean energy dissipation rate for the suggested 
methods does not depend much on the step size. Figure [5] shows the trajec- 
tory of the pendulum in the x-y-plane. Notice that that Heun's method does 
not "stabilize" around a limit cycle, which is the case for Algorithms |3.1f|3.3[ 
Instead, the trajectory drifts outwards in an exponentially increasing fashion. 
This is because Heun's method does not preserve a nearby relative equilibrium, 
whereas the analyzed methods do. We now continue with a discussion of relative 
equilibrium for the elastic pendulum. 

Relative equilibria For each momentum value /i e g = TidS*^ = M there is 
a corresponding relative equilibrium for the elastic pendulum, corresponding to 
rotation with constant angular velocity about the vertical axis. 

Turning to cylindrical coordinates {p, (j), z) on =S, and corresponding co- 



tangent lifted momenta {pp,p^,Pz) on T*J3, the governing equations (I6I take 
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Energy diagrams for elastic pendulum 




20 40 , 60 80 100 



Figure 4: Energy diagrams for the elastic pendulum, using Algorithms |3.1H3.3| 
and Heun's method, with e — 0.01 (upper) and e = 0.1 (lower). The energy 
oscillates for Algorithms |3. lH3?3l but shows that qualitatively correct behavior, 
independently of the step size. For Heun's method, the result depends heavily 
on the step size, and is qualitatively incorrect, since energy is increasing. 
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Projected 2: y trajectory for elastic pendulum 

3S,2S,RKS {h = 0.1, e = 0.1) Heun {h = 0.1, e = 0.1) 

2- 




-2-10 1 2 

3S,2S,RKS {h = 0.1, e = 0.01) 





-2-10 1 2 

3S,2S,RKS {h = 0.4, e = 0.01) 




Figure 5: Trajectories in the x-y-plane for the elastic pendulum, using Algo- 
rithms |3]T}|3^ and Heun's method. With Heun's method, the trajectory does 
not stay bounded, but drifts outwards. With Algorithms |3.1fj3.3| the solution 
approaches a limit cycle (corresponding to a relative equilibrium). The exact 
correct limit cycle is not attained, but instead a modified one, in accordance 
with the analysis in Section [6] 
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Modified relative equilibria for Algorithms |3.1|43.3| 



4^ 




Figure 6: Near preservation of relative equilibria for Algorithms |3.1f|3.3| 
with e = 0.1 (the results between the methods are inseparable). Initial data 
is given on the exact relative equilibrium cycle. The trajectory drifts towards a 
modified relative equilibrium cycle 0(/i^) away from the exact cycle. The result 
verifies the analysis in Section [6] 



the form 
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ep 
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Notice that Jz{q,p) = p^,- Due to the 5^-invariance of g,d,V, these equa- 
tions reduce (in accordance with co-tangent reduction theory) to a system on 
T{^/S^). We may use the coordinates {p, z,Pp,Pz), in which case the reduced 
equations are obtained by simply removing the equations for (f) and p^, and 
replacing with the momentum value p. 

The equilibrium on the reduced space, corresponding to momentum value p, 
is given for initial data 7(0) = {p, z,pp,pz) that fulfill the relations 



kp{ 



-1 - 



kz 



1 



mp- 

- mg 



3=0, Pp = 
= 0, p,^ 0. 



Thus, the corresponding relative equilibrium solution 7(<) = cxp(t^) • 7(0) is 
given by 

Pp{t) = 

tp 



'pit) 


= P(0), 


< m 


= 0(0) 




= z(0). 



to(p(0))2' 



p4t) 

Pz{t) 



P 

0. 
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Since the action is given by exp(t^) • {p,(j),z) = {p,<j) + t^,z), it holds that 
^ = ij,/{mp{0)^) for the equilibrium solution. 

The x-y-trajectory using Algorithms |3. 1|-|3.3| and e — 0.1, with initial data 



Qo = (0, 1.0392304845413263, -0.6), Po = (1.3416407864998736, 0, 0), 

corresponding to exact relative equilibrium for momentum value fi = p^, is given 
in Figure |6] Again, the result between the methods is inseparable. Notice that 
the trajectories slightly drift off the true relative equilibrium, toward a modified 
relative equilibrium. The rank of d is one, meaning that the dissipation is de- 
generate even after the reduction (it is zero if (z, p) and {pz^Pp) are orthogonal). 
Thus, we cannot apply Proposition |6.1 to explain the stability of the modified 
relative equilibrium for Algorithm |3.2f 3.3[ However, notice that the reduced 



Hamiltonian is strictly convex in a neighborhood of the critical point [7^(0)]. 
Thus, by Theorem |6.2| the behavior seen in Figure |4] is fully explained for Al- 
gorithm |3.1[ 



7 Conclusions 

For Hamiltonian problems on standard form perturbed by Rayleigh damping, 
three numerical integration schemes were analyzed. When the dissipation pa- 
rameter is zero, all the methods reduce to the Stormer-Verlet scheme. Based on 
backward error analysis, we gave result on: 



Asymptotically correct dissipation rate of the modified Hamiltonian (The- 
orem 



4.3 and Theorem 4.6 1 



• Monotone decay of the modified Hamiltonian for small enough step sizes 



(Theorem 4.6 1 



Preservation of symmetries, and conservation of momentum (Proposi- 
tion 



5.2) 



Near preservation of relative equilibria (Theorem 6.2 1 



The theoretical results were verified by numerical examples of a damped planar 
pendulum, and a damped elastic 3-D pendulum. 
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